;+
; NAME:
;       MAKE_MW_PLOT
;
; PURPOSE:
;       Plot Galactocentic positions of sources atop the pretty
;       Spitzer artist's rendering.
;
; CATEGORY:
;       Plotting
;
; CALLING SEQUENCE:
;       MAKE_MW_PLOT
;
; INPUTS:
;       NONE
;
; OPTIONAL INPUTS:
;       NONE
;
; KEYWORD PARAMETERS:
;       NONE
;
; OUTPUTS:
;       NONE
;
; OPTIONAL OUTPUTS:
;       NONE
;
; COMMON BLOCKS:
;       NONE
;
; MODIFICATION HISTORY:
;
;       Created:  11/21/12, TPEB -- Initial version.
;       Modified: 12/20/12, TPEB -- Turned into a PRO and added
;                                   (simple) documentation for filing
;                                   in the repository.
;       Modified: 03/11/13, TPEB -- Moved Bar circle to 4.0 kpc a la
;                                   Benjamin:2005.
;       Modified: 03/18/13, TPEB -- Added B&W version for published
;                                   paper to avoid paying $$$.
;
;-

PRO MAKE_MW_PLOT
  
     ;; Add to the plot our stuff...
     mw = omni_read_conffile('./conffiles/galactic_params.conf')
     dpdfs = omni_read_conffile('./conffiles/dpdf_params.conf')
     d = dindgen(dpdfs.NBINS)*dpdfs.BINSIZE + dpdfs.BINSTART
     R0 = MW.R0
     
     restore,'./local/BGPS_V2_pvec.sav',/ver
     ;; restore,'./local/BGPS_V2_emaf.sav',/ver
     ;; n = n_elements(emaf)
     ;; eind = where(constrain.emaf)
     ;; pvec = pvec[eind]
     ;; constrain = constrain[eind]
     
     pvec.stat.duse /= 1.d3
     
     ;; Restrict consideration to "well-constrained" sources...
     cind      = where(constrain.post,n)
     ;; emaf      = emaf[cind]
     pvec      = pvec[cind]
     constrain = constrain[cind]
     
  bw=[0,1]
  suf=['','_bw']
  
  ;; Run for color & BW
  FOR colt=0,0 DO BEGIN
     
     
  ;; Make the plot for both the annotated image and the blank image
  outfn = ['annot','blank']
  infn  = './emaf_paper/MW_Spitzer_'+['A','B']+suf[colt]+'.png'
  
  FOR jj=0,1 DO BEGIN
     
     myps,'./grs_paper/MW_plot_'+outfn[jj]+suf[colt]+'.eps',BW=bw[colt],$
          xsize=3.1,ysize=3.1,thick=1.5
     
     
     ;; the Spitzer Images have 134.379 pix/kpc in the 5600x5600 version
     ;; We're using the 1500x1500 version:
     pixper = 134.379 * 1500./5600.
     ;; Also, cropping to 33/56*1500 = 884 pixels  (Down to 680 pix now)
     ;; crop = 3300./5600.*1500. * 680./884.
     crop = 3300./5600.*1500. * 884./884.
     width = crop / pixper
     
     ;; Include a fudge-factor for scaling to properly place sun at 
     ;;   8.5 kpc, edge = width/2.
     edge = (width/2.)*1.065 *(r0/8.5d3) ; Now place at 8.3 kpc
     range = [-edge,edge]
     
     ;; !x.thick=3
     ;; !y.thick=3

     cs = 0.55
     
     ;; Plot the image using cgImage
     cgImage,image,filename=infn[jj],/keep_aspect,charsize=cs,$
             xtit='x!dgal!n [kpc]',xrange=range,yrange=range,$
             ytit='y!dgal!n [kpc]',color='Black',/interpolate,$
             background='White',/axes,position=[0.11,0.07,0.99,0.99]
     
     
     ;; ;; l=30 line (through end of the Bar)
     rr0 = r0/1.d3
     elr = (findgen(11)+1.)*15.*!const.dtor
     slopes = -1.d / tan(elr)
     sunx = 0.
     suny = -rr0
     xoff = [0.5,0.8,0.5,2.5,2.5,0.5,0,0,0,0.9,0]
     yoff = [-1.0,-1.0,-0.2,-1.2,-0.4,0.3,0,0.2,0,0.2,0]
     FOR ll=0,n_elements(elr)-1 DO BEGIN
        
        eld = elr[ll]*!const.rtod
        
        IF (eld GT 90) && (fix(eld) MOD 30 NE 0) THEN CONTINUE
        print,eld
        
        xp = eld GE 90. ? sunx : -2.d * rr0 * cos(elr[ll]) * sin(elr[ll])
        yp = slopes[ll]*xp + suny
        
        x2 = !x.crange[0]
        y2 = slopes[ll]*x2 + suny
        
        IF y2 GT !y.crange[1] THEN BEGIN
           y2 = !y.crange[1]
           x2 = (y2-suny) / slopes[ll]
        ENDIF
        
        IF y2 LT !y.crange[0] THEN BEGIN
           y2 = !y.crange[0]
           x2 = (y2-suny) / slopes[ll]
        ENDIF
        
        cgPlots,[xp,x2],[yp,y2],$
                linestyle=0,thick=3,color='bisque'
        
        cgText,x2+xoff[ll],y2+yoff[ll],charsize=0.7*cs,'l = '+$
               string(eld,format="(I0)")+cgSymbol('deg'),$
               orientation=eld-90
     ENDFOR
     
     ;; Solar Circle and Tangent Circles
     ells = findgen(1001)/1.e3*180.
     scl = findgen(1001)/1.e3*360.
     xells = -(R0/1.d3 * cos(ells*!dtor)) *  sin(ells * !dtor)
     yells = + (R0/1.d3 * cos(ells*!dtor)) * cos(ells * !dtor) - R0/1.d3
     xsc   = -R0/1.d3 * cos(scl * !dtor)
     ysc   = -R0/1.d3 * sin(scl * !dtor)
     xbar  = -4.0 * cos(scl * !dtor)
     ybar  = -4.0 * sin(scl * !dtor)
     cgOplot,xells,yells,linestyle=1,color='Black',thick=2
     cgOplot,xsc,ysc,linestyle=1,color='Black',thick=2
     cgOplot,xbar,ybar,linestyle=2,color='Blk7',thick=2
     
     ;; Make x,y position for plotting
     x = -pvec.stat.duse[0] * sin(pvec.glon*!dtor)
     y = -R0/1.d3 + pvec.stat.duse[0] * cos(pvec.glon*!dtor)
     
     ;; Plot all the BGPS Objects!!!1!
     cgPlots,x,y,psym=16,symsize=0.2,color='Black'
     
     ;; Sun and Galactic Center
     cgPlots,0,0,psym=7,color='BLK7',symsize=0.6,thick=3
     cgPlots,0,-R0/1.d3,psym='filledcircle',color='bisque',symsize=0.6
     cgPlots,0,-R0/1.d3,psym='opencircle',color='black',symsize=0.6,thick=1
     cgPlots,0,-R0/1.d3,psym='filledcircle',color='black',symsize=0.15,thick=1
     cgText,0.1,0.4,'GC',charsize=cs*0.5,color='Black'
     cgText,0.1,-8.0,'Sun',charsize=cs*0.5,color='Black'
     
     ;; Plot BeSSeL trigonometric Parallaxes?
     restore,'./ancillary/bessel_parallaxes.sav',/ver
     d_px = 1.d / bessel.px
     omni_lbd2rz,bessel.l,bessel.b,d_px*1.d3,brg,bzg,XGAL=bxg,YGAL=byg
     bxg /= 1.d3
     byg /= 1.d3
     print,m4_stat(bxg,byg)
     FOR qq=0,n_elements(bessel)-1 DO BEGIN
        IF byg[qq] LT !y.crange[0] THEN CONTINUE
        cgPlots,bxg[qq],byg[qq],psym='filledsquare',symsize=0.25,color='yellow'
        cgPlots,bxg[qq],byg[qq],psym='opensquare',symsize=0.25,color='black',$
                thick=0.25
     ENDFOR
     
     ;; ;; Add HRDS?
     ;; restore,'./ancillary/HRDS_catalog.sav',/ver
     ;; omni_lbd2rz,hrds.glon,hrds.glat,hrds.dsun*1.d3,hrg,hzg,XGAL=hxg,YGAL=hyg
     ;; hxg /= 1.d3
     ;; hyg /= 1.d3
     ;; FOR qq=0,n_elements(hrds)-1 DO BEGIN
     ;;    IF hyg[qq] LT !y.crange[0] THEN CONTINUE
     ;;    cgPlots,hxg[qq],hyg[qq],psym='filleddiamond',symsize=0.25,color='magenta'
     ;;    cgPlots,hxg[qq],hyg[qq],psym='opendiamond',symsize=0.25,color='black',$
     ;;            thick=0.25
     ;; ENDFOR
     
     
     ;; Label spiral arms, as desired...
     cgText,-6.2,-1.5,'Sgr Arm',charsize=cs*0.7,$
            alignment=0.5,orientation=90.,color='black'
     cgText,-5.7,3.9,'Perseus Arm',charsize=cs*0.7,$
            alignment=0.5,orientation=40.,color='Black'
     cgText,-10.5,-1.0,'Outer Arm',charsize=cs*0.7,$
            alignment=0.5,orientation=80.,color='Black'
     cgText,-0.7,-2.5,'End of Bar?',charsize=cs*0.7,$
            alignment=0.5,orientation=0.,color='Black'
     cgText,+1.8,-5.7,'Sct-Cen Arm',charsize=cs*0.7,$
            alignment=0.5,orientation=0.,color='Black'
     
     
     cgAxis,xaxis=0,/xst,xtickformat='blank_axis'
     cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
     cgAxis,yaxis=0,/yst,ytickformat='blank_axis'
     
     
     
     myps,/done
     
  ENDFOR
  
ENDFOR
END
